# 谷歌瓦片下载

import re
import math
import random
from pathlib import Path

import urllib.request as ur
from threading import Thread, Lock

MAP_URLS = {
    "google": "http://mt{server}.google.com/vt/lyrs={style}&hl=zh-CN{offset}&src=app&x={x}&y={y}&z={z}",
    "amap": "http://wprd02.is.autonavi.com/appmaptile?style={style}&x={x}&y={y}&z={z}",
    "tencent_s": "http://p3.map.gtimg.com/sateTiles/{z}/{fx}/{fy}/{x}_{y}.jpg",
    "tencent_m": "http://rt0.map.gtimg.com/tile?z={z}&x={x}&y={y}&styleid=3"
}

COUNT = 0
mutex = Lock()


# ------------------wgs84与web墨卡托互转-------------------------

# WGS-84经纬度转Web墨卡托


def wgs_to_mercator(x, y):
    y = 85.0511287798 if y > 85.0511287798 else y
    y = -85.0511287798 if y < -85.0511287798 else y

    x2 = x * 20037508.34 / 180
    y2 = math.log(math.tan((90 + y) * math.pi / 360)) / (math.pi / 180)
    y2 = y2 * 20037508.34 / 180
    return x2, y2


# Web墨卡托转经纬度


def mercator_to_wgs(x, y):
    x2 = x / 20037508.34 * 180
    y2 = y / 20037508.34 * 180
    y2 = 180 / math.pi * \
         (2 * math.atan(math.exp(y2 * math.pi / 180)) - math.pi / 2)
    return x2, y2


# -------------------------------------------------------------


# ---------------------瓦片地址到墨卡托---------------------------------
'''
东经为正，西经为负。北纬为正，南纬为负
lon经度 lat纬度 z缩放比例[0-22] ,对于卫星图并不能取到最大，测试值是20最大，再大会返回404.
山区卫星图可取的z更小，不同地图来源设置不同。
'''


# 根据WGS-84 的经纬度获取谷歌地图中的瓦片坐标


def wgs84_to_tile(lon, lat, z):
    '''
    gps              tile
         ^              -------->
         |              |
    -----|---->  =>>    |
         |              v

    Get google-style tile cooridinate from geographical coordinate
    lon : Longittude
    lat : Latitude
    z : zoom
    '''

    def isnum(x):
        return isinstance(x, int) or isinstance(x, float)

    if not (isnum(lon) and isnum(lat)):
        raise TypeError("lon and lat must be int or float!")

    if not isinstance(z, int) or z < 0 or z > 22:
        raise TypeError("z must be int and between 0 to 22.")

    if lon < 0:
        lon = 180 + lon
    else:
        lon += 180
    lon /= 360  # make lon to (0,1)

    lat = 85.0511287798 if lat > 85.0511287798 else lat
    lat = -85.0511287798 if lat < -85.0511287798 else lat
    lat = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180)
    lat /= 180  # make lat to (-1,1)
    lat = 1 - (lat + 1) / 2  # make lat to (0,1) and left top is 0-point

    num = 2 ** z
    x = math.floor(lon * num)
    y = math.floor(lat * num)
    return x, y


def tile_to_mercator(tile_x, tile_y, z):
    length = 20037508.3427892
    sum = 2 ** z
    mercator_x = tile_x / sum * length * 2 - length
    mercator_y = -(tile_y / sum * length * 2) + length
    return mercator_x, mercator_y


# -----------------------------------------------------------


# ---------------------度分秒转换------------------------------

def dms2dd(degrees, minutes, seconds, direction):
    dd = float(degrees) + float(minutes) / 60 + float(seconds) / (60 * 60)
    if direction == 'E' or direction == 'N':
        dd *= -1
    return dd


def dd2dms(deg):
    d = int(deg)
    md = abs(deg - d) * 60
    m = int(md)
    sd = (md - m) * 60
    return [d, m, sd]


def parse_dms(dms):
    parts = re.split('[^\d\w]+', dms)
    lat = dms2dd(parts[0], parts[1], parts[2], parts[3])
    return (lat)


# -------------------栅格像素到墨卡托----------------------------

def imagexy2geo(trans, u, v):
    '''
    根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
    :param trans: GDAL地理转换矩阵
    :param v: 像素的行号
    :param u: 像素的列号
    :return: 行列号(v, u)对应的投影坐标或地理坐标(x, y)
    '''
    px = trans[0] + u * trans[1] + v * trans[2]
    py = trans[3] + u * trans[4] + v * trans[5]
    return px, py

# -----------------GCJ02到WGS84的纠偏与互转---------------------------


def transformLat(x, y):
    ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * \
          y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
    ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 *
            math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(y * math.pi) + 40.0 *
            math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 *
            math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
    return ret


def transformLon(x, y):
    ret = 300.0 + x + 2.0 * y + 0.1 * x * x + \
          0.1 * x * y + 0.1 * math.sqrt(abs(x))
    ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 *
            math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(x * math.pi) + 40.0 *
            math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 *
            math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
    return ret


def delta(lat, lon):
    """
    Krasovsky 1940
    //
    // a = 6378245.0, 1/f = 298.3
    // b = a * (1 - f)
    // ee = (a^2 - b^2) / a^2;
    """
    a = 6378245.0  # a: 卫星椭球坐标投影到平面地图坐标系的投影因子。
    ee = 0.00669342162296594323  # ee: 椭球的偏心率。
    dLat = transformLat(lon - 105.0, lat - 35.0)
    dLon = transformLon(lon - 105.0, lat - 35.0)
    radLat = lat / 180.0 * math.pi
    magic = math.sin(radLat)
    magic = 1 - ee * magic * magic
    sqrtMagic = math.sqrt(magic)
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * math.pi)
    dLon = (dLon * 180.0) / (a / sqrtMagic * math.cos(radLat) * math.pi)
    return {'lat': dLat, 'lon': dLon}


def outOfChina(lat, lon):
    if (lon < 72.004 or lon > 137.8347):
        return True
    if (lat < 0.8293 or lat > 55.8271):
        return True
    return False


def gcj_to_wgs(gcjLon, gcjLat):
    if outOfChina(gcjLat, gcjLon):
        return gcjLon, gcjLat
    d = delta(gcjLat, gcjLon)
    return gcjLon - d["lon"], gcjLat - d["lat"]


def wgs_to_gcj(wgsLon, wgsLat):
    if outOfChina(wgsLat, wgsLon):
        return wgsLon, wgsLat
    d = delta(wgsLat, wgsLon)
    return wgsLon + d["lon"], wgsLat + d["lat"]


# --------------------------------------------------------------


# ---------------------下载器相关-----------------------------
agents = [
    'Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/60.0.3112.101 Safari/537.36',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/532.5 (KHTML, like Gecko) Chrome/4.0.249.0 Safari/532.5',
    'Mozilla/5.0 (Windows; U; Windows NT 5.2; en-US) AppleWebKit/532.9 (KHTML, like Gecko) Chrome/5.0.310.0 Safari/532.9',
    'Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US) AppleWebKit/534.7 (KHTML, like Gecko) Chrome/7.0.514.0 Safari/534.7',
    'Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/9.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/10.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.20 (KHTML, like Gecko) Chrome/11.0.672.2 Safari/534.20',
    'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/534.27 (KHTML, like Gecko) Chrome/12.0.712.0 Safari/534.27',
    'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.1 (KHTML, like Gecko) Chrome/13.0.782.24 Safari/535.1'
]


class Downloader(Thread):
    # multiple threads downloader
    def __init__(self, index, count, urls, datas, update):
        # index 表示第几个线程，count 表示线程的总数，urls 代表需要下载url列表，datas代表要返回的数据列表。
        # update 表示每下载一个成功就进行的回调函数。
        super().__init__()
        self.urls = urls
        self.datas = datas
        self.index = index
        self.count = count
        self.update = update
        self.pos = None

    def download(self, url):
        HEADERS = {'User-Agent': random.choice(agents)}
        # print(url)

        err = 0
        while err < 4:
            full_url = url.format(server=err)  # 如果下载失败，则切换服务器
            header = ur.Request(full_url, headers=HEADERS)
            try:
                # print('download ', full_url)
                data = ur.urlopen(header).read()
            except:
                err += 1
                print("Bad network link.", full_url)
            else:
                [x, y, z] = self.pos
                path = "titles/{}/{}/".format(z, y)
                Path(path).mkdir(exist_ok=True, parents=True)
                path += "{}.jpg".format(x)
                save_img(path, data)
                return self.pos, data
        # raise Exception("Bad network link.")
        print("Bad network link.", url)
        return None

    def run(self):
        for i, (x, y, z, url) in enumerate(self.urls):
            if i % self.count != self.index:
                continue
            self.pos = [x, y, z]
            self.download(url)
            # self.datas[i] = self.download(url)
            if mutex.acquire():
                self.update()
                mutex.release()


def downTiles(urls, multi=10):
    def makeupdate(s):
        def up():
            global COUNT
            COUNT += 1
            print("\b" * 45, end='')
            print("DownLoading ... [{0}/{1}]".format(COUNT, s), end='')

        return up

    url_len = len(urls)
    datas = [None] * url_len
    if multi < 1 or multi > 20 or not isinstance(multi, int):
        raise Exception(
            "multi of Downloader shuold be int and between 1 to 20.")
    tasks = [Downloader(i, multi, urls, datas, makeupdate(url_len))
             for i in range(multi)]
    for i in tasks:
        i.start()
    for i in tasks:
        i.join()

    return datas


def save_img(path, bit):
    with open(path, "wb") as f:  # wb是写二进制
        f.write(bit)


def geturl(source, x, y, z, style, offset):
    """
    Get the picture's url for download.
    style:
        m for map
        s for satellite
    source:
        google or amap or tencent
    x y:
        google-style tile coordinate system
    z:
        zoom
    """
    if source == 'google':
        offset = '&gl=CN' if offset else ''
        furl = MAP_URLS["google"].format(
            server='{server}', x=x, y=y, z=z, style=style, offset=offset)
    elif source == 'amap':
        # for amap 6 is satellite and 7 is map.
        style = 6 if style == 's' else 7
        furl = MAP_URLS["amap"].format(x=x, y=y, z=z, style=style)
    elif source == 'tencent':
        y = 2 ** z - 1 - y
        if style == 's':
            furl = MAP_URLS["tencent_s"].format(
                x=x, y=y, z=z, fx=math.floor(x / 16), fy=math.floor(y / 16))
        else:
            furl = MAP_URLS["tencent_m"].format(x=x, y=y, z=z)
    else:
        raise Exception("Unknown Map Source ! ")

    return furl


def getTilesByBBox(bbox, zoom):
    """
    bbox依次输入左上角的经度、纬度，右下角的经度、纬度，缩放级别，地图源，输出文件，影像类型（默认为卫星图）
    获取区域内的瓦片并自动拼合图像。返回四个角的瓦片坐标
    """
    x1, y1, x2, y2 = bbox
    pos1x, pos1y = wgs84_to_tile(x1, y1, zoom)
    pos2x, pos2y = wgs84_to_tile(x2, y2, zoom)
    lenx = pos2x - pos1x + 1
    leny = pos2y - pos1y + 1
    print("Total number：{x} X {y}".format(x=lenx, y=leny))

    tiles = []
    for y in range(pos1y, pos2y + 1):
        for x in range(pos1x, pos2x + 1):
            tiles.append([x, y])

    return tiles, [pos1x, pos1y, pos2x, pos2y], [leny * 256, lenx * 256]


def getUrlsByTiles(tiles, zoom, source='google', style='s', offset=False):
    urls = []
    for x, y in tiles:
        url = geturl(source, x, y, zoom, style, offset)
        urls.append([x, y, zoom, url])
    return urls


def downloadTitle(zoom, bbox):
    # TODO 设置不同服务商的地图
    source = 'google'
    # TODO 设置不同类型 s为影像
    style = 's'
    # zoom = 18
    offset = True
    # bbox = [120.09198099000002, 27.176422685000115,
    #             120.09819027500012, 27.173431772500086]
    tiles, tile_bbox, image_shape = getTilesByBBox(bbox, zoom)
    urls = getUrlsByTiles(tiles, zoom, source, style, offset)
    downTiles(urls)
    print("\nDownload Finished！ Pics (w,h) is (%d,%d) Mergeing......" %
          (image_shape[0], image_shape[1]))


def downloadTitles():
    # TODO 等级范围
    zoom = [18, 18]
    # TODO bbox范围
    bbox = [120.09198099000002, 27.176422685000115,
            120.09819027500012, 27.173431772500086]
    for z in range(zoom[0], zoom[1] + 1):
        downloadTitle(z, bbox)


if __name__ == "__main__":
    downloadTitles()
